In this document ...
In [1]:
# import header
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import climlab
from climlab.model.ebm import EBM_seasonal
from climlab.solar.orbital_cycles import OrbitalCycles
from climlab.surface.albedo import StepFunctionAlbedo
from climlab.domain.field import global_mean
An EBM model instance is created through
In [2]:
# model creation
ebm_model = EBM_seasonal(A=204, B=2.17, D=0.79, a0=0.3, a2 = 0.01, ai=0.6, Tf=-2)
The model is set up with an albedo feedback subprocess:
In [3]:
# print model states and suprocesses
print ebm_model
In [4]:
ebm_model.integrate_converge()
global_mean(ebm_model.timeave['Ts'])
Out[4]:
In [5]:
ebm_model.diagnostics['icelat']
Out[5]:
In [6]:
ebm_model.time
Out[6]:
The model integration through a different set of orbital parameters is done by a function called OrbitalCycles from the climlab.solar.orbital_cycles module.
We want to integrate the model for orbital parameters from the last 5,000 years. To save computation time we speed up the process by the factor 10.
In [7]:
# run for 5,000 orbital years, but only 500 model years
experiment = OrbitalCycles(ebm_model, kyear_start=-23, kyear_stop=-10, segment_length_years=10., orbital_year_factor=100.)
In [8]:
ebm_model.time
Out[8]:
In [9]:
from climlab.solar.orbital import OrbitalTable
In [32]:
kyears = np.arange(-50,0,0.1)
table = OrbitalTable()
orb = table.lookup_parameters(kyears)
In [42]:
fig = plt.figure( figsize = (9,9) )
ax1 = fig.add_subplot(311)
ax1.plot(kyears, orb['ecc'])
ax1.axvline(-23,color='r')
ax1.axvline(-10,color='r')
ax1.set_title('Eccentricity $e$', fontsize=18 )
ax2 = fig.add_subplot(312)
ax2.plot(kyears, orb['obliquity'])
ax2.axvline(-23,color='r')
ax2.axvline(-10,color='r')
ax2.set_title('Obliquity (axial tilt) $\Phi$', fontsize=18 )
ax3 = fig.add_subplot(313)
ax3.plot(kyears, orb['ecc'] * np.sin( np.deg2rad( orb['long_peri'])))
ax3.axvline(-23,color='r')
ax3.axvline(-10,color='r')
ax3.set_title('Precessional parameter $e \sin(\Lambda)$', fontsize=18 )
ax3.set_xlabel( 'kyears before present', fontsize=14 )
Out[42]:
In [12]:
experiment.T_segments_annual[45,:]
#experiment.orb_kyear
#experiment.kyear_start
#experiment.segment_length_years * experiment.orbital_year_factor /1000.
Out[12]:
In [13]:
#import matplotlib.pyplot as plt
#from mpl_toolkits.mplot3d import Axes3D
# creating plot figure
fig = plt.figure(figsize=(15,10))
lat = ebm_model.lat
start = experiment.kyear_start
stop = experiment.kyear_stop
interval = experiment.segment_length_years * experiment.orbital_year_factor / 1000.
time = np.arange(start, stop, interval)
# Temperature plot
ax1 = fig.add_subplot(111)
ax1.contour( time, lat, experiment.T_segments_annual, cmap='RdBu')
#ax1.set_xticks([-90,-60,-30,0,30,60,90])
#ax1.set_xlim([-90,90])
#ax1.set_xlabel('latitude')
#ax1.set_ylabel('surface temperature (degC)')
#ax1.grid()
plt.show()